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CN ■ Abstract 

l__3 ■ In this paper, we propose an efficient acquisition scheme for GPS receivers. It is shown that GPS 

signals can be effectively sampled and detected using a bank of randomized correlators with much 
fewer chip-matched filters than those used in existing GPS signal acquisition algorithms. The 
latter use correlations with all possible shifted replicas of the satellite-specific C/A code and an 
exhaustive search for peaking signals over the delay-Doppler space. Our scheme is based on the 
recently proposed analog compressed sensing framework, and consists of a multichannel sampling 
structure with far fewer correlators. 

The compressive multichannel sampler outputs are linear combinations of a vector whose sup- 
port tends to be sparse; by detecting its support one can identify the strongest satellite signals 
in the field of view and pinpoint the correct code-phase and Doppler shifts for finer resolution 
during tracking. The analysis in this paper demonstrates that GPS signals can be detected and 
acquired via the proposed structure at a lower cost in terms of number of correlations that need 
to be computed in the coarse acquisition phase, which in current GPS technology scales like the 
product of the number of all possible delays and Doppler shifts. In contrast, the required number 
of correlators in our compressive multichannel scheme scales as the number of satellites in the field 
of view of the device times the logarithm of number of delay-Doppler bins explored, as is typical 
for compressed sensing methods. 
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1. Introduction 

Nowadays, Global Positioning System (GPS) chips are ubiquitous, and continue to be embedded 
in a variety of devices. A GPS device allows to determine its location with about 3 meters accuracy, 
by measuring the propagation delay of signals transmitted by the set of GPS satellites in the 
field of view (FOV) of any receiver located on the surface of the earth, which typically requires 
measurements from at least four satellites [l|. 

Conventionally, the signal that arrives at the receiver is downconverted, match-filtered and 
oversampled at a fast rate. Subsequently, the receiver acquires enough (at least four) strong signals 
by exploiting the orthogonality of the distinct coarse/acquisition (C/A) codes used in GPS signaling 
at each satellite [2j]. However, due to the unknown propagation delays, the samples obtained are 
misaligned in time and frequency and therefore, it is vital to pinpoint the code-phase in order to 
decode the navigation data correctly [2][3j and use the time-delay information for pseudo-range 
computation. Furthermore, each of the satellites contributes a component of the received GPS 
signal that is characterized by a distinct Doppler offset [4|], due to the unequal relative velocity of 
satellite and receiver, as well as the offset of the different local oscillators at the GPS receivers. 
In general, time-frequency synchronization as well as signal detection is tackled in GPS receivers 
during the acquisition/ detection stage via a parallel search over the binned delay-Doppler space 
across all the satellite C/A codes [5] [6]. 

In many practical scenarios, signals might arrive at the receiver with multipath components 
instead of the line of sight (LOS) component [3] [6] [7]. Constructive and destructive superposition 
of randomly delayed and faded replicas, leads to distorted correlation peaks. This is usually tackled 
in the tracking stage [2|] that follows the acquisition/detection stage, by using an early-late receiver. 
Such a receiver compares the energy of a symbol period in the first half from the early gate to the 
energy in the last half from the late gate so that the receiver can synchronize the signals accordingly. 
Furthermore, many approaches, in addition to the early-late structure, have been proposed to better 
mitigate the effects brought by multipath, including (but not limited to) the Narrow Correlator 
[g], Multipath Eliminating Technique (MET) [9j, and Multipath Estimating Delay Lock Loop 
(MEDLL) [10l | . These methods differ in their capabilities to remove multipath errors, specifically 
at low signal-to-noise ratio (SNR) and/or in the presence of interference. In this work, we consider 
the general signal model that considers multipath effects and propose an acquisition scheme that 
coarsely captures significant paths for each active satellite, with its corresponding code-phase and 
Doppler. The tracking stage that further resolves the estimates of delay-Doppler pairs as well as 



the multipath components is beyond the scope of this paper. 

As described above, the acquisition and detection of GPS signals is usually performed sequen- 
tially. First, the strongest signals coming from the satellites are detected by searching a binned 
delay-Doppler space via exhaustive correlations that pinpoint the correct coarse timing information 
and frequency offsets. After acquisition and detection, the signal is locked and the device enters 
the tracking stage that tackles fine synchronization and multipath error mitigation in order to 
despread, demodulate and decode the navigation data correctly in real-time. However, this acqui- 
sition/detection scheme can be computationally intensive due to the large number of correlations, 
and especially the exhaustive search for peaks over the binned delay-Doppler space across all the 
satellite signals with distinct C/A codes. For example, the maximum Doppler shift in a GPS signal 
is typically within [— 10kHz, 10kHz] and the search step size is usually 500Hz while the maximum 
delay can run up to a C/A code length 1023. In this case, the 2-D delay-Doppler peak is found 
by comparing the outputs of 1023 x 41 ~ 4 x 10 4 correlators for each satellite, which is a heavy 
computation burden. 

Paper contributions: In order to scale down the operations and hardware requirements, we 
propose a simple and efficient acquisition scheme based on the recently developed compressed 
sensing (CS) framework [131 ] and its extension to analog signals 151 ] . The multichannel samplers in 



[151 ] are constructed as a randomized linear combination of the duals of all the generators, where the 
generators in this case correspond to the satellite-specfic C/A code waveforms. In our context, we 
show that the the duals of the generators are well approximated by the generators themselves. This 
alleviates one of the most difficult aspects in the practical application of [15], namely, the physical 
implementation of the dual filters, by exploiting properties of the spread spectrum sequences that 
are in the GPS standard. Thanks to this interpretation, the proposed multichannel samplers can 
be viewed as performing independent random projections of all correlators outputs. The resulting 
set of compressive measurements are then used together to recover the peaks located sparsely over 
the delay-Doppler space, which is a jointly sparse recovery problem with infinite input vectors and 
infinite measurement vectors (IMV). The continuous-to-finite (CTF) method introduced in [18] 
effectively reduces the IMV problem to a finite multiple measurement vector (MMV) system with 
jointly sparse inputs, which can be solved efficiently using the Reduce MMV and Boost (ReMBo) 
technique proposed in [18[ , or other MMV approaches [l6(] [13] . 

The paper is organized as follows. Section 2 describes the general model for GPS signals. 
Section 3 re-interprets existing GPS acquisition schemes from a sampling point of view. In order 



to scale down the computations and hardware requirements, Section 4 introduces the analog CS 
framework. In Section 5 we further reduce the general solution to a set of simple compressive 
samplers by utilizing the structure of GPS signals. Numerical results are shown in Section 6 to 
demonstrate the effectiveness of our proposed acquisition scheme, followed by a complexity analysis 
given in Section 7. Finally the paper is concluded in Section 8. 

2. GPS Signal Model 

The signal transmitted by the satellites is a direct sequence spread spectrum (DS-SS) signal 
modulated onto LI and L2 frequencies at 1575.42MHz and 1227.60MHz respectively. In commercial 
GPS systems publicly available to civil users, the DS-SS signal received at the user end is carried 
on LI frequency from all the available launched satellites. Equivalently, the baseband signal from 
the ith satellite is transmitted as 

Si {t) = J2 d i[n}Mt-nT), i = l,---,I (1) 

where (j>i{t) is a spreading waveform determined by a satellite-specific spreading code and {dj[ n ]}nez 
is the navigation data sent by the ith satellite with a symbol period of T, containing its time stamp, 
orbit location and relevant information entailed for positioning the receiver. 

More specifically, the waveform (fti(t) is determined by the ith. satellite's C/A code {sj[ m ]} as 

A/-1 
&(*)= ^2si[m]g(j:-mT c ), i = !,-■■, I (2) 

m=0 

where g(t) is a wideband short pulse. For simplicity, we assume that g(t) has a flat spectrum of 
bandwidth $l g = 2ttL/T c (typically L = 1) approximated with error e g (cu) 

G(u) = [1 + e g (u))] rect 27r L/T c (w), (3) 

where e g (co) specifies the deviation from the flat spectrum with_| ||e 5 (w)|| <C 1. Due to the periodicity 
of the C/A code, T = MT C . The C/A code {sj[?^]} is a pseudo-random binary sequence of length 
M that contains A^ maximum length sequence (MLS) or Gold sequence of length Mq = 1023 



lr The notation || ■ || refers to the L2 norm of a function ||e 9 (w)|| — \//^° |e9(^)| 2 dtJ 



transmitted with a chip period T c = 977.5ns, which implies M = NMq. In fact, by the GPS 
transmission standards we have N = 20 for the GPS Ll-C/A signal, i.e. T = 20ms. 

The correlation properties of the spreading code are vital in the recovery of spread spectrum 
signals. Denote the cross-correlation between different C/A code as 

M-l 

Ri'iH = jj Yl Si '^ m ~ w l s iN- ( 4 ) 

m=0 

When M is large, the Gold sequences or MLS sequences are orthogonal between different satellites 
and approximately orthogonal between different shifts |7j. This is indicated by the flat and 2ir/T c - 
periodic cross spectral density 

M-l 

<M e " Tc )- E Ri'M<r' 1 ™ Tc = stf - %] + e iV H, (5) 

u=-M+l 

where the error function e^ t j(u;) is also 27r/T c -periodic witro ^^^(o;)]! <C 1. This flat property plays 
an essential role in simplifying the design presented later in this paper. 



After downconversion, the signal at the receiver can be modeled as 

/ R 
x ( f ) = Yl Y, h hMt - n, r )e^^ + v(t), (6) 



where {/ii, r }r=i, •■-,/? are the multipath channel taps with delays {Ti r }r=l,— ,R and Doppler shifts 
{uji i r} r =i l — ,r from the ith satellite to the receiver, and v(t) is the Additive White Gaussian Noise 
(AWGN) with variance a 2 . Combined with the signal model ([1]), the signal x(t) is represented as 

/ R 

X(t) = YjYjJ2 °i,rM&(* ~nT- T^)e lU] ^ + v (fi t 
n£Z i=l r=l 

where aj ir [n] = hi tr di[n]. In the coarse acquisition phase, it is typically assumed that the delays are 
integer multiples of the chip duration n r = qi )T T c with qi tr £ Q and the Doppler shifts are integer 
multiples of the frequency search step u>i r = ki r Aui with hi r £ )C, where the sets Q and /C define 



The norm here is defined as ||ei',i(w)|| — J_/t l e »',*( w )l ^ due to the periodicity. 



the delay-Doppler space. This leads to the following discretized signal model 

/ R 



x(t) = Y,Y,J2 a i,rHMt -nT- q itr T c )e ik '- A ^ + v(t). (7) 

neZ i=l r=l 

3. Standard GPS Acquisition Scheme 

The main task of the acquisition stage is to detect the correct code-phase q = {<&,?• }j=i'... V 
and Doppler shift k = {fe,r}i=l'- j across the delay-Doppler space and recover the sequence 
{ai i? .[n]}j = i r .. /, among which the strongest set I of satellites (\I\ > 4) are picked for the purpose of 



triangulation 
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[4|. Note that the sequence {aj )r [n]} nG g includes the attenuation of the channels 
between the satellites and the receiver. Therefore, its magnitude indicates the strength of the signal 
received and only the strong ones are acquired by the receiver. In general, the magnitudes of those 
acquired i £ I are significantly greater than those i ^ X, making the coefficients a, jr [n] sparse due 
to the wide difference in signal strength. 

3.1. Exhaustive Search via Matched Filtering (MF) 

Conventionally, the acquisition and detection of strong satellite signals is achieved by correlating 
the incoming signal x(t) with a bank of match-filters (/>«(£) 's that are separately modulated by 
carriers {e *}fcgx: and shifted in time {4>i(t — qT c )} g€ Q. In this way, the paths corresponding to 
peaks in the magnitude of aj j? .[n] can be found in the delay-Doppler binned-space Q x fC for each 
satellite corresponding to its C/A code. 

This approach can be viewed as sampling with a set of filters, followed by uniform sampling at 
time t = nT, as depicted in Fig. 1. The sampling kernels of this equivalent structure are given by 
4>i,k,q{t) = 'Pii't — qT c )e lkAujt , for all i = 1, ■ ■ ■ , I, k G K. and q G Q. The sampled output in each 
channel is equal to 

Zi,k,q[n] = (x(t), <t>i,k,q(t - nT)). (8) 

In the Fourier domain, we have 

**. (e'- T ) - k t, •*», (» - 2 4) x (" - 2 4) ■ < 9 > 

eez ^ ' ^ ' 

where $*fc„(w) and X{u) are the Fourier transforms of <f>i,k,q(—t) and x(t) respectively. Note 
that the summation over £ G 7L in ([9]) depends on the bandwidth of the filter (j)i t k,q(—t)i where as 
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Figure 1: Exhaustive Matched Filtering (MF) Approach 

mentioned in Section 2 the bandwidth of g(t) is Vt g = 2ttLM/T. Therefore, the summation becomes 
finite from £ = to £ = LM — 1 over to € [—tt/T, tt/T]. From ([7]), we can express X(lj) as 



/ R 



X(u>) = YjJ2 Ai ' r ( e ^ T ) ^ - ^,rAw)e- i(w - <: '- rAu, ^ rTc + V {e iujT ) , 



(10) 



=1 r=l 



where we defined A i<r (e iwT ) = YLn&L a^ r \n]e~' in{ ^- k ^^^ T . Substituting (JTUJ) into ©, and denot- 
ing by z (e lajT ) the length-/|/C||Q| column vector whose (i,k,q)th element is Zi t k,q (e lwT ), and by 
en (e lajT ) the length- R column vector of {A ir (e lwT )} r=1 „ for the ith data stream, we can write 



z (e iwT ) = M^( W) k, g)a (e iwT ) + v (e iwT ) 



(11) 



over the domain u S [—tt/T, tt/T]. The derivation is identical to the development in [lj] and is 
therefore omitted. Here a (e lwT ) = [a^ (e lwT ) , • • • , a^ (e lwT )]^ is a length-/i? vector containing 
the DTFT of all the data sequences {ai,r[ n ]}nez an d Mw,j,(w, k, q) is an i"|/C||Q| x IR matrix with 
[(i, k,q), (i,r)]th element 



LM-l 



r™ i , m 1 ^ \* ( 27r A 



$ 



^i^i.r ?Qi,r I 



2vrA 



(12) 



The component v (e lwT ) = [•••, ^j,fc,g (e ia;T ) , • • • ] T is the filtered noise by matched filters (genera- 
tors) {4>i,k,q{p)} i=\ '. q . i and therefore has a cross-spectral density matrix H vv (e lajT ) = a 2 'M. ( j )( j ) (ui , /C, Q), 
where M^(w, /C, Q) is the Gram matrix of all the generators defined by 



LM-l 



[M^( W) /c, Q)]^.^.^ = ^ E $ ^w ( - - ^ ) *'■*•" 1 - — ' 



>(i' ,k',q'),(i,k,q) r n 



£=0 



2tt£ 



2tt£\ 



(13) 



Exploiting the specific choice of sampling kernels and the structure of M00(w,fc, q) and M<^(u;,/C, Q), 
we can further analyze the output samples z (e 1 ) as stated below. 

Theorem 1. Suppose that the following conditions hold, 

CI) the pulse shaping filter has a spectrum G{uj) = [1 + e(cj)] rect 27r £/T c ( w ) with error e g (oj); 

C2) the C/A code cross spectral density is SVj (e lwTc ) = 5[i' — i] + ei> : i(u>) with error e^^u); 

C3) the frequency search step size is chosen as Aoj = 2tt j/T and j € Z + . 

If the error functions satisfy \\e g (uj)\\ <C 1 and ||ej' j(u;)|| <C 1 for any i',i = 1, • • • ,/, then the Gram 
matrix of all the generators {4>i{t — qT c )e t } i=1 ] q . j satisfies 

M # ( W ,/C,Q) = LM + E( W ), (14) 

where E(w) is bounded perturbation matrix satisfying \\\^{^)](i' k' q') (i k q)\\ = ®(X) ^ LM and 
the filtered noise samples have a cross-spectral density matrix H vv (e lwT ) = a 2 [LMI + E(u;)]. The 
output samples z[n] = [■ ■ ■ , Zikq[n], ■ ■ ■] at each of the kernels 4>i t k,q(t) ccn be written as 



Zi,k,q[n\ 



LMa iyr [n} + 0(1) + v ijkiq [n], q = q^ r and k = k iyT 
0(1) + Vi t k,q[n], otherwise, 



where Vi t k,q[n] is the time-domain filtered noise sample and 0(1) <C LM is some bounded perturba- 
tion error with LM being the processing gain on the signal-to-noise ratio. 

Proof See Appendix A. □ 

Note that the frequency step size Auj = 2nj/T corroborates the fact that for standard commer- 
cial GPS systems, the step size is usually 2tt x 500 rads/s which fits the analysis here by choosing 
j = 10. Also, we can see that the output Zi^k,q[n] at each sampler represents the correlation be- 
tween the matched filters and the incoming signal, which is proportional to the magnitude of a^j-fn] 
and corrupted by noise. Assuming large enough processing gain LM and small enough noise, the 
delay-Doppler pairs {t^ t = qi, r T c }lZi '... j an d {wi, r = fc^rAw}^'.'"'; can be found by the location 
of the peaks/dominant entries in Zj fc g[ro]- The strongest set of satellite signals can then be detected 
by comparing the values in Zik q [n] so that a subset X of the satellite signals are locked and passed 
onto the tracking stage for finer extraction. If we ignore the noise for a moment, then 2j fc o[w] is 
sparse in the sense that for each value n it contains only a small number of non-zero entries. 



3.2. Compressive Multichannel Acquisition 

Although effective, this conventional approach taken by standard GPS receivers performs ex- 
haustive correlations (MF approach) that requires abundant samples from a large number of cor- 
relators I|Q||/C|. This task can be computationally expensive and demanding on the hardware and 
memory resources. Assuming a maximum channel delay spread of r max = QT C and Doppler shift 
of l w max| = KAuj, the total number of correlators is 2IQK. For example, the maximum Doppler 
shift is typically ±10kHz. Assuming a delay spread up to code length r max = MT C , then with a 
frequency grid of 500Hz, the total number of correlators needed becomes 24 x 1023 x 42 ~ 10 6 . 

Therefore, it is highly desirable to scale down the computational complexity and power con- 
sumption of a user GPS device by performing less correlations while sustaining its capability to 
pinpoint the signal timing and Doppler information during acquisition. By observing the correlation 
outputs in the vector z[n] = [• • • , 2» fc «|Vi]) ■ ■ ■ ] T , it can be seen that only few of the dominant entries 
are useful. Our goal is to exploit the underlying sparsity in the signal model to design an acquisition 
scheme that requires far fewer correlators. Instead of tackling the problem from a match-filtering 
viewpoint as in standard GPS, we look at the problem from an analog CS perspective [la], which 
is one of the main contributions of this paper. 

The analog CS design outlined in [15| requires a small number of samplers (only twice the 
sparsity 2|X|i? in a noiseless setting), and hence gives rise to substantial practical savings as analyzed 
later in Section 7. However, the solution [15] is given in the frequency domain and in general does 
not admit a tractable form in time domain, which makes it hard to implement in practice. Another 
contribution of this work lies in further exploiting the structure of GPS signals so that the sampling 
kernels are easy to implement. The outputs from the compressive samplers can then be used to 
solve the sparse recovery problem of locating the dominant /peak values reflected in the vector z[n], 
for example, using the method in [181 ] . 

Before we go into the details of our design, we start by describing the analog CS framework 
[151 ] . In Section 5 we further develop and simplify the general solution to fit our problem. 

4. Compressed Sensing of Analog Signals 

The exhaustive MF approach in standard GPS receivers acquires delays and Dopplers by di- 
rectly exposing the sparse structure in the output samples z[n] obtained from a large number of 
correlators. In order to reduce the number of correlators while retaining the ability to correctly 
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Figure 2: Compressive Multichannel Sampling (CS) 

identify the peaks of z[n], it is possible to directly measure a compressed version of z[n] at the 
samplers outputs and recover that sparse structure instead, by employing analog CS techniques. 

4-1. General Model for Analog Compressed Sensing (CS) 

The signal model in ([7]) does not reflect any sparse structure, since it is expressed by a set 
of deterministic generators </>j(i)'s defined by unknown parameters qi jT and ki tr . The sparsity we 
exploit is the sparsity of delay-Doppler pairs pinpointed by the peaks/dominant entries in z[n] over 
the entire delay-Doppler space Q x /C for each user i = 1, ■ ■ ■ ,1 that is informative in acquiring the 
signal. Using a dictionary {4>i(t — qT c )e lk t }k£K,q£Q, the signal can be equivalently expressed by 

x{t ) = E E E E y^MUt -nT- qT c )e ikA ^ + v(t), 
nez i=l keKqeQ 



where 



Vi,k,q[n] = < 



&j,r[/i], Q = qi,r and k = ki t ., 
0, otherwise. 



(15) 



Note that the sparsity of y[n] = [■ • ■ , yi t k,q [n], •••] is identical to that of z[n] in the noiseless setting. 
Indeed, for each i = 1,- ■ ■ ,1 there are altogether R dominant coefficients {yi,k,q[n]}keK,q£Q that 
correspond to the original coefficients {ai jr [n]} r= i t ... t R and select the correct code-phase qi tT and 
Doppler shifts fcj jT .. Let the support of y[n] be S, then the support S contains the code-phase and 
Doppler information for acquisition, with a sparsity of |<S| = \I\R. The aim of analog CS is to 
exploit this sparsity in acquiring x{t) using fewer correlators. 
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4-2. General Solution of Compressive Samplers 



As the scheme of [15[ uses a set of compressive samplers ip p (—t), p = 1, • • • ,P<€. I\IC\\Q\ to 
obtain minimal measurements, from which the sparse vector y[n] can be recovered. As depicted in 
Fig. 2, the samples at the output of ip p {—t) at t = nT are given by 

c p [n]±{x(t),^ p (t-nT)). (16) 

Similar to the mathematical manipulations in Section 3, the system equation can be re- written as 

c (e iwT ) = M^(w, K, Q)y (e iwT ) + w (e^ T ) , (17) 

where M^(u>, /C, Q) is a P x J|/C| | Q| matrix with similar structure to (|12p and the notation w (e 1 ) 
is used to distinguish the noise component from the previous method in standard GPS. It has been 
proven in [15|] that in a noiseless setting, simply twice the sparsity P = 2\I\R is needed for successful 
recovery of the sparse vector y[n], if il> p (— t)'s are chosen properly. For noisy scenarios, the necessary 
number of channels P is larger than the minimum, and evaluated numerically; in any case, it is 
much smaller than that required by the standard scheme, as we will demonstrate in Section 6. 

This reduction is obtained by appropriately choosing a set of randomized correlators \l/(u;) = 
[ x I / i(a;), • • • , ^p(lo)] t . A general expression of the compressive samplers is given in [15[ as 

*( w ) = BM:!( w ,/c,Q)$( w ,x;,g), (is) 



where B is a sensing matrix satisfying certain coherence properties [13J (e.g., Gaussian random 
matrix or partial DFT matrix [13], or an appropriate deterministic binary matrix [19(), and 
*(w,/C, Q) = [••• ,$ i (w-A;Au;)e" i ( CJ " feA ^' ?Tc ,---] T isalength-/|/C||Q| vector containing the Fourier 
transforms of the generators {<fii(t — qT c )e *} i= i ] q . j ■ With this choice of *S?(u)), it can be shown 
that M^(w,/C, Q) = B. Since B is independent of frequency u), transforming (|T7|) into the time 
domain, the samples can be written as 

c[n] = By[n] + w[n], n G Z. (19) 

The vectors {y[w]} are jointly sparse since they all share the same sparsity pattern. To find y[n], we 
can convert (|19p to a finite MMV problem using the continuous-to-finite (CTF) technique developed 

11 



in [18(]. Specifically, we first find a basis for the range space of {c[n]} by computing the covariance 
matrix R cc and decomposing it as R cc = CC . Here C can be chosen as the eigenvectors of R cc 
multiplied by the square-root of the corresponding eigenvalues. Then, the support of y[n], n £ Z 
can be obtained by solving C = BY, where Y is the sparsest matrix satisfying the measurement 
equation. This problem can be treated using various MMV sparse recovery techniques [16|][l7t]. In 
our simulations, we use the ReMBo algorithm developed in [18J . Finally the support of y[n] is 
obtained by taking the union of the supports of the columns in the matrix Y. Once the support of 
y[n] is recovered, the acquisition of correct delay-Doppler pair is automatically achieved by locating 
the dominants/peaks in the vector y[n] of (I15j) . 

Remark: As verified in Section 6, as the number of correlators P increases, the acquisition 
performance improves significantly. Solving the MMV problem requires collecting multiple mea- 
surement vectors {c[n]}, while the standard GPS scheme can either employ information for a single 
measurement z [n] in (J8J) or further leverage the processing gain over multiple measurements {z [n] } . 
For the proposed compressive acquisition scheme, if a single vector measurement is used to recover 
the sparse vector y[n] using greedy methods or ^i-norm based methods, the performance will de- 
grade as shown in Fig. [3] but not significantly. Therefore, there is a trade-off between the number 
of observations c[n], the number of acquisition channels P as well as the accuracy of the acquisition 
in comparison with the standard GPS scheme. 

5. Simplified Randomized Correlators 

The method proposed in [15] depends on the ability of physically implementing the sampling 
kernels in (|18j) . Therefore, we explore the structure of the matrix M^^w, fC,Q) to provide practical 
insights on the design of such filters. 

Corrolary 1. Suppose that the conditions (C1)-(C3) and the requirement on the error functions 
in Theorem 1 hold. Then the sampling kernels can then be chosen as the randomized correlators 

I 

Mt) = E E E W,«)&(* - qT c )e ikA ^, p = 1, ■ ■ ■ , P. (20) 

i=i kefCqeQ 

Proof. From (|18f) we have the general solution of the compressive samplers 

*(u)=BM£(u,K,Q)*(u,K,Q). (21) 
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Algorithm 1 Compressive Multichannel Acquisition 



(1) Construct P compressive sampling kernels as 



Mt) =EEE b P,(^ q )Mt - qT c )e ikA "\ (23) 

j=i keicqeo. 



where [B] p ^ j. „ = b pi k q is a sensing matrix that satisfies certain coherence properties such as RIP 13 [. 

(2) Apply the set of compressive sampling kernels \&(u;) and arrive at measurements 

c[n] = By[n] + w[n], neZ. 

(3) Solve the jointly sparse recovery problem as in [18( to recover the support of y[n]. 

(4) Once the support of y[n] is available, the delay-Doppler pairs are determined by the support q = q,^ r 
and k — ki r as in (11511 . 



According to the result in Theorem 1, using Taylor expansion on the matrix inverse M , , (w,/C, Q) 
and ignoring high order terms scaled by 1/LM <1, we can approximate the inverse by 

1 + Tm e ^) = i ~ m E ^ + <w EV) " <W eV) ■ ■ • M '' (22> 

where the last approximation comes from the fact that E(u>) contains negligible elements. Therefore, 
the compressive samplers can be chosen directly as 9{u) = B<&(w,/C, Q), which leads to the time- 
domain expression in the corrolary. □ 



The The filter responses of (|20|) can be precomputed, and these P channel outputs are sampled 
every T = MT C to produce the test statistics that are going to be used in lieu of the coefficients 
z[n] in Theorem 1. 

Remark: Note that although the samples are taken at 1/T, the physical implementation of 
the compressive multichannel filtering operation is likely to require digital processing at the chip 
rate 1/T C . Nevertheless, it is possible that a wise choice of the coefficients of the matrix B can 
further help reduce computations while maintaining the identifiability of the parameters. Analysis 
of this approach goes beyond our current scope. What we can certainly claim is that the number 
of computations is now controlled by the parameter P, rather than by the number of possible 
generators that span all possible delays Q and Dopplers K. In fact, the sampling Kernels are 
precomputed and used online. This is likely to reduce cost of computation, access to memory and 
storage. The performance of the compressive multichannel sensing structure degrades gracefully as 
P decreases, giving designers degrees of freedom to choose a desirable operating point. 
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6. Numerical Results 

In this section, we run numerical simulations to demonstrate the proposed CS acquisition scheme 
in GPS receivers. In the simulation, \X\ = 4 out of I = 24 satellites asynchronously transmit 
C/A signals that are received by the GPS devices, where the codes are length-M = NMq Gold 
sequences with N = 20 and Mq = 1023. A total of n = 50 navigation data bits are sent at the 
rate of \/T = 50Hz (i.e., T = 20ms). The transmit filter is modeled by a finite length pulse shaping 
filter g(t) = \/T~ c sm.c{t/T c ) when |t| < T g , and g(t) = otherwise. The length T g is sufficiently 
large so that the response of the pulse in the frequency domain remains approximately flat, i.e. 
G(u) «rect 27r/Tc (w). 

To reduce the simulation overhead without incurring a loss of generality, we assume that our 
statistical model for the channel consists of uniformly distributed delays, Tj jr ~ U(0, r max ) that are 
bounded by a maximum delay spread of r max = 20T C ; and of Doppler shifts, \u}i >r \ < w mM that 
are uniformly distributed, Wj jr . ~ U{— w max ,u> max ) over a frequency range delimited by u> max /27r = 
2.5kHz. The channel gains are hi r ~ CJ\f(0, 1), with a multi-path propagation having R = 2 
paths per satellite. In order to identify fractional delays with a half-chip accuracy, the functions 
0i,fc,g(O = <fri(t ~ qT c )e * are chosen with a half-chip spacing q = 0, 1/2,1,- •• such that the 
resolution of At = T c /2 is achieved, and with a frequency resolution of Aw = 10 x 2tt/T that 
corresponds to steps around 500Hz when T = 20ms. It follows that \Q\ = [r max /Ar] + 1 = 41 
and \JC\ = 2 [w max /Au;] + 1 = 11. For simulation purpose, the sensing matrix B is generated as 
a random binary matrix (while in practice it can be chosen as a deterministic binary matrix to 
simplify the implementation of correlators [19||). 

In all simulations, the attenuated components with distinct delays from each of the satellites 
are acquired by a number of P = {80, 120, 240, 360, 480} channels, in contrast to the traditional 
24x41xll~lx 10 4 . The performance is illustrated in terms of success rate and average Root 
Mean Square Error (RMSE), respectively, in Fig. [3] and Fig. [4j The success rate of acquisition is 
the probability P(X = I) of the proposed scheme to determine the strongest \I\ = 4 signals, which 
is shown in the figure against the number of channels P and the SNR. The conditional RMSE is an 
average error between the true delay-frequency parameters and those associated to the strongest 
paths of the correctly identified satellites 



RMSE average (g) 



A 



\ |{xnx}| ^ ^ ' 



14 



Satellite Identification Success 




Satellite Identification Success Rate 



50 
100 














































X,", 


















250 
300 

350 













0; 

t 



-40 -35 -30 -25 -20 
SNR(dB) 



Figure 3: Satellite identification rate P(Z = I) using a 
single measurement c[l] for a CS receiver v.s. the MF 
receiver for P = {120,240,360,480} (above), and for 
P = {10, 20,..., 400} (below) 
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Figure 4: Satellite identification rate P(Z = I) using 
multiple measurements {c[l], c[2], . . . , c[50]} for CS re- 
ceiver v.s. the MF receiver for P = {80, 120, 240, 360} 
(above), and for P = {10, 20, . . . , 400} (below) 



where (Tj,u;i) = (ri, r * , U)i >r * ) with r* = argmax r£ r 1 m. |^i,r| 2 and 



(%,/ci) = argmax|zj ifci q[n]| = argmax|y iifci(? [n] 
qeQ,k€tc qeQ,keic 



(24) 



are the delay-frequency index pairs of strongest path associated to the ith satellite. Similarly, the 
average RMSE for the Doppler is 



RMSE. 



average 



(fc) 



N 



1 



\{ini}\ 



^ (/CjAw - LUi) C < 



ie{ini} 



Although the compressive acquisition scheme suffers from a compression loss, both FigJ3]and FigU] 
highlight its ability to perform closely as the traditional MF. When P > 80 and SNR > —25 dB 
the active satellites I can be identified satisfactorily which leads to great savings (less than 1% of 
the original 1 x 10 4 ). 

The figures above illustrate acquisition performances using a single set of measurements c[0] 
against that using multiple sets of measurements {c[n]}^'L 1 . Using a single measurement suffers 
from a performance loss (—10 dB for P = 120 at the rate of approximately 0.8). In fact, by reducing 
n, the accuracy of z[n] and consequently the sensitivity, degrade. Furthermore, it can be seen from 
Fig. [3] that the required number of channels P has to be raised to 480 (less than 5% of the original 
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Delay RMSE Performance 



Delay RMSE Performance 




Figure 5: Delay estimation (above) and Doppler esti- 
mation (below) performance of the CS, with n = 1 and 
P = {120, 240, 360, 480} compared against the the MF 
receiver 




Figure 6: Delay estimation (above) and Doppler (esti- 
mation) performance of the CMA, with multiple mea- 
sures (n = 50) and P = {80,120,240,360}, compared 
against the the MF receiver also processing n — 50 
measures 



1 x 10 4 ) to achieve a reliable rate that approaches the MF result. 

A similar trend is also visible on the conditional RMSE curves for both single (Fig. [5]) and 
multiple (Fig. E) modes (-12 dB for P = 120 when RMSE(g) « 2 and RMSE(fc) « 2 • 1(T 3 ). At 
high SNR the performance is limited by the presence of a systematic error due to the modeling 
mismatch from the quantized parameters. At low SNR, instead, the error is bounded by the length 
of the search interval QT C . Once again the CS method closely approaches the MF performance, 
especially when n = 1. 



7. Complexity Analysis 

The complexity of the acquisition algorithm is due to two aspects: 1) storage requirement and 
2) computational complexity. We provide here a brief analysis of the complexity of the proposed 
CS scheme against traditional MF scheme. To make a fair and practical comparison, we assume 
that the implementation is done in the digital domain and we use the Lk eme i-tap digitized version 
of the sampling kernels {Vv (*)}»= l (and also {<l>i,k,q(t)}i=i '. q . j for the traditional case). 
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7.1. Storage and Processing Requirement 

The difference in storage results from two sources, one is the storage for the digital kernel taps 
and the other is the outputs of the sampling kernels used for peak recovery, both of which are 
proportional to the number of sampling kernels. Furthermore, the processing overhead per unit of 
time for these stored values scales proportionally with the storage requirement as well. 





sampling kernels 


output samples 


CS Receiver 


■* x -^kernel 


O(P) 


MF Receiver 


I|/C||Q x L kernc i 


0(I\K\\Q\) 



It is clear that the proposed compressive acquisition scheme handles less data, which facilitates the 
pipelining of the algorithm and also relieves the burden of storage. 



7.2. Computational Complexity 

The difference in computations stems from the correlations and the search for the peak. The 
number of operations in performing correlations is proportional to the number of sampling kernels, 
while the peak recovery is different for the two approaches, depending on how the sparse recovery 
(proposed CS structure) and the exhaustive search (MF structure) are implemented. Here we 
further compare the two architectures by their number of operations that are necessary to identify 
the delay-Doppler pairs (I24p . In this practical analysis, the compressed samples c[n] are obtained 
by post-processing of the digitally sampled versions of x(t ) at the chip rate and processed using a 



greedy algorithm Orthogonal Matching Pursuit (OMP) 2l|. Note that using analog implementation 
in the acquisition can further bring down the complexity in terms of processing. 

We introduce a vector x[n] of M dimensions, whose mth entry is {x[n]} m = x(nT + mT c ), to 
digitally capture and compress one instance of the signal according to 

c p [n] = (x[m],ip p [m — nM\). (25) 

For the MF receiver, instead, we assume an oversampling ratio of 2 to achieve half chip accuracy, i.e., 
At = T c /2, and downsize the filterbank array. The sequence x[n] is partitioned into 2 sub-sequences 
{xi[n],X2[n]}, of M samples each, whose mth element is {xj[n]} m = x(nT + mT c + (i — l)T c /2), 
i = 1,2. A typical filter model would process the stream of 2M samples sequentially, however to 
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CS Receiver 



Complexity 



Remarks 



Digital compression c p [n] 
Covariance R cc (optional 3 ) 
SVD of R cc (optional 3 ) 
Residual update 
Inner products 
Maximum projection 
Least-Squares projection 
Stopping criterion 



0(nMP) 

0(nP 2 ) 

0(n 2 P) 

0(n\S\ 2 ) 

0{nPI\K\\Q\\S\) 

0(|5|log(/|/C||Q|)) 

o(\s\ 3 ) 

0{nP\S\) 



Eq. ([25]) 

MMV mode 

n<P [22], MMV mode 

(0MP.1) 

(0MP.2) 

(OMP.3) 

(OMP.4) 

(OMP.5) 



MF Receiver 



Complexity 



Remarks 



Correlations Zi^, q [n\ 
Path selection 
Accumulation 



0{nMI\K\\Q\) 
0(nIR log(|/C||Q|)) 
0(nI\lC\\Q\) 



Eq. 



MMV mode 



Table 1: Complexity breakdown for the proposed CS and traditional MF acquisition using n sets of measurements 

emulate the block processing nature of the CS receiver and avoid CPU cycles that would further 
delay the execution of the algorithm, we let the 2 sub-sequences be processed concurrently. 

All the arithmetic operations, starting from x[n], necessary to detect the \I\R vector elements 
are recorded and listed in Table [TJ The table outlines both single and multiple (MMV) modes 
for both the MF and CS schemes, and a breakdown of the OMP recovery algorithm adopted 
by the CS receiver. This popular algorithm seeks the S (with \S\ = \I\R) non-zero elements of 
the sparse vector y[n] by sequentially choosing dictionary elements that better correlate with the 
observations c[n]. At every iteration the current estimate is subtracted from the observation vector 
(OMP.l) and the residual projected onto the dictionary elements (OMP.2). Then, the dictionary 
element linked to the largest coefficient (OMP.3) is retained and removed from the dictionary. 
The updated set of coefficients is obtained by projecting c[n] onto the subspace formed by the 
set of atoms that were removed from the dictionary (OMP.4). The algorithm stops when either a 
maximum number of iterations is reached or when the norm of the residual falls beyond a predefined 
threshold (OMP.5). 

Path selection refers to identifying the support of a certain vector for pinpointing the active 
components (delay-Doppler pairs). For both CS and MF, it is tightly coupled to the sorting 
algorithm being implemented and therefore, we only list its average computational complexity 
rather than the number of comparators. 

When the representation of y[n] is sufficiently sparse, i.e. for GPS applications |<S| <C 7|/C||Q|, 
the number of operations needed to digitally compress x[n] into c[n] and to project the residual 
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Average Synchronization Runtime (single) 




Figure 7: Average runtime to evaluate {I, q, k} from a received observation vector y[n] for the compressive 
scheme, with n = 1 (above) and n = 50 (below) as a function of P = {20, 40, . . . , 480}, and compared against 
the MF receiver. Each curve was run separately on a 64-bit i7 920 CPU running at 2.67 GHz. 



of each OMP iteration onto the dictionary (OMP.2) dominate the overall complexity of the CS 
receiver, leading to an order of 0(nP(M + 7|/C||Q||<S|)). On the other hand, the number of oper- 
ations for the MF are mainly determined by the number of additions to compute the correlations, 
leading to C(nM/|/C||Q|). For hardware implementation this is attractive since the filterbank pro- 
cessing does not require complex multiplications. However, a similar saving can be added to the 
CS receiver by appropriately designing B such that {V'pb 71 ]} a ls° nas ±1 elements. 

The comparison between the dominant terms results in a CS to MF complexity ratio {P/I\K\ \ Q\ + 
P\S\/M) that favors the former and emphasizes the complexity savings. In fact, one should in fact 
expect P|«S| <C M and P <C /|/C||Q|, which shows that the CS gains by removing its dependency 
on the length of the C/A sequence. This trend is also highlighted in Fig. [7] by the average CPU 
time spent while executing the steps described in Table [TJ 

When n > 1 the ratio remains unchanged since all the additional steps (Table [1]) for the ReMBo 
technique require marginal increase of operations. When compared to n = 1, the MF spends more 
CPU time to accumulate the correlation outputs whereas the CS receiver experiences a reverse 
trend. The additional efforlo spent to evaluate R cc is compensated by less operations within the 
OMP algorithm, and results in a gain in efficiency as highlighted in Table [TJ 

In general, knowing a priori the order |<S| the CS receiver has an advantage over the MF, which 



3 Note that in practice, the covariance and SVD computation can be optional by directly choosing a set of mea- 
surements {c[n]} and solve the MMV instead. 
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is true and practical in GPS sytems because the order of number of active satellites in the field 
of view is actually known. However, the MF approach always explores and ranks all the |/C||Q| 
dimensions for every satellite before selecting \X\R. 

8. Conclusions 

We proposed a compressive multichannel acquisition scheme for GPS receivers. The reduction 
is achieved by choosing randomized linear combinations of all the MFs, which leads to great savings 
in practice. As shown in the analysis and numerical results, our scheme can efficiently recover the 
unknown delay-Doppler pairs using significantly fewer correlators than those needed in a standard 
GPS receiver. Regardless of the sparse recovery algorithm, the acquisition performance improves 
gracefully with the increase of acquisition channels and the number of observations. Therefore, 
although the proposed scheme has a performance loss in terms of RMSE and success rate compared 
to the standard GPS scheme, it provides a design tool to trade-off complexity and performance 
gracefully that can be useful to scale down the cost and energy consumption of GPS chips. 
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Appendix A. Proof of Theorem 1 



In this proof, we prove the structure of the matrices M<^(w, k, q) and M^(w, /C, Q), which will 
lead to the results of the proposed theorem. 

Let <f>i, k , q (-t) = cpi{t-qT c )e lkAwt . Denote by 9 itk>q (u) = $t(w - kAco) e - i( - u - kAu ^ c the Fourier 
transform of <pi t k,q(—t) ■ Using 3>i(cu) = G(oj) X^ m =o Si[m]e~ imulTc together with the spectrum 
G{uj) = [1 + e g (oj)] rect 27r L/T c (w) and ignoring higher order perturbations C(|e 9 (cj)| 2 ), we can write 
the [(«', k, q), (i, r)]th entry of the matrix M^(w, fc, q) of (|12p over w G [— 7r/T, 7r/T] as 



[M H (uj,k,q)] (i , k 



1 



LM-l 



(i',k,q),(i,r) ~ T / / 

A/-1 M"-l 



(A.l) 



J^ J^ 4[m']s i [m]e iw(m '- m)Tc e- i t( m '- m )e- iA ^( m ' i - m ^) + 0(e 9 ( w )). 



rn'=0 m=0 



With a change of variable u = m — m', we can re- write the double summations over m and m! as 



M-l M-l 



Yl Yl 4[™> ! He iu(m '" m)Tc e- i ^ (m '- m) e" iA * (mVm ^ r) 



m'=0 m=0 
M-l 



M-l 



J2 M-^^s^m- u} Si [m}e- iAuT ^ k - k ^ m e --(-^-¥)^, 



(A.2) 



u=-M+l 



m=0 



^R^^Ujk—ki^ 



where Ri'i[u,k — k^ r ] is a pseudo-correlation between the sequence {sj/[m]} and {sj[m]} being 
perturbed by phase-shifts determined by the mismatch of the Doppler shift e - lAulT c( k - k i,r) m _ Based 
on the results in [20] and and taking into account that |sj[m]| = 1, it can be shown that 



±\jjf j ti, /C fXij 



M-l 

— Y s i'[ m - u]si{m], k = fej,, 



M 



m=0 
M-l 



M ^ 
m=0 

0(1/M), 



U j=- Q,k / fej, r , 



(A.3) 



where 0(1/M) is some small perturbation. In particular, it is desirable for the pseudo-correlation 
Ri'i[u, k — fej r ] to decay rapidly over k such that dominant values only appear when there is a 
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frequency component kAoj = ki ;r Aco. By choosing Aw = 2irj/T, j £ Z + , we have 

M-l M-l 

}_ J2 e -^T c{k - Kr)m = _(_ £ e -i^(*-*.,r) = 6 [k - k hr ], (A.4) 

m=0 m=0 

which results in a simplified expression of the pseudo-correlation as follows 

R n [u, k - k i>r ] = Rn[u]8[k - k i>r ] + 0{l/M), 

where Ri'i[u] is the C/A code cross-correlation in (jl|). Therefore, the matrix entry contains signifi- 
cant values only if k = ki^ r , 

LM-l 



f|\/r ,,(,, V. n \] — _Ju(q-qi,r)T c -ik i:r (q-q lir )AujT c ST^ p -i^(q-qi,r)T c 

[ w H<t>K U ' K ' < i)\(i,ki ir ,q) ) {i,r)— T e e / . 

1=0 

M-l 

x Y, R^[u\e-' m ^- kA ^^) Tc +0{e g {u)) + 0{l). 



■^ 



ui-kAui-^i )T C 



Furthermore, using the spectrum Sm (e 1 c ) = S[i' — i] + £j',i(w) and ignoring higher order pertur- 
bations O(\ei' t i(oj)\ 2 ), the non-zero entries of the matrix M^(a;, k,q) are explicitly written as 

LM-l 

e=o 
+ O{e g {oj)) + O(e i , )i {oj))+O{l). 

With T = MT C , we use the property 

LM-l 
1 £ e -i¥(«-«.r)Tc =Wb _ ftfr] (A . 5) 

to further express the non-zero entries of rVLw,(w, fc, q) at i' = i, k = ki >r and g = qj jr 

[M^(w,fc, </)],. fe , )( . r) =LM + C»(e 9 (a;)) + 0(e M (w)) + C(l), w e [-7r/T )7 r/T]. (A.6) 
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On the other hand, the matrix M^,(w,/C, Q) can be expressed element-wise as 



LM-l 



Similarly, the expression is significant only when k = k! 

LM-l 

[M^u, /C, Q)]^^^) = pU-M'Ptor*)***. £ e i^(g- 9 ')T c + 0(eg(w)) + 0(et , >i(w)) + (l). 
When g = g', according to ()A.5p . we have 

[M^( W ,/C, G)] (i , )fcjg))(i)M) = LM + 0(e») + 0(e iV (u;)) + 0(1). 
Since the error functions satisfy ||cp(w)|| <C 1 and ||€i',t(w)[| <C 1, the results in Theorem 1 follow. 
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